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Abstract 

In this article, the so-called "Nystrom method" is tested to compute optimal quantizers of Gaus- 
sian processes. In particular, we derive the optimal quantization of the fractional Brownian motion 
by approximating the first terms of its Karhunen-Loeve decomposition. 

A numerical test of the "functional stratification" variance reduction algorithm is performed with 
the fractional Brownian motion. 
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Introduction 



Let (ft, A, P) be probability space, and E a reflexive Banach space. The norm on E is denoted | • |. 

The quantization of a random variable X, taking its values in E consists in its approximation by 
a random variable Y taking finitely many values. The resulting error of this discretization is the LP 
norm of \X — Y\. Minimizing this error, with a fixed maximum cardinal of Y(£Y) yields the following 
minimization problem: 

min{||X- Y|| p ,y : fi -> E measurable, card(Y"(fi)) < N} . (1) 

This problem, was first considered for signal transmission and compression issues. More recently, quan- 
tization has been introduced in numerical probability, to devise quadrature methods [16J, solving multi- 
dimensional stochastic control problems [S] and for variance reduction [3]. Since the 2000's, the infinite 
dimensional setting has been investigated from both theoretical an numerical viewpoint, especially in the 
quadratic case [12]. One elementary property of a L 2 optimal quantizer is the stationarity: E[X|Y] = Y. 
If X is a bi-measurable stochastic process on [0, T] verifying J T E[|X t | 2 ](it < oo, it can be considered 
random variable valued in the Hilbert space H = L 2 ([0,T]). In [12], it is shown that in the centered 
Gaussian case, linear subspaces U of H spanned by TV-stationary quantizers correspond to principal 
components of X, in other words, are spanned by eigenvectors of the covariance operator of X. Thus, 
the quantization consists first in exploiting its Karhunen-Loeve decomposition (e^,A^) >1 < 

If d (N) is the dimension of the subspace of L 2 ([0,T]) spanned by Y(SY), the quantization error 
Cn(Y) writes 



(X) 2 = J2 Af + ejy (g) AT(0, Af ) for m > d N (X). (2) 



J > 771 + 1 




eN( 



e N {X) 2 < X f + eN I A f ) ] for 1 < "* < d N {X). (3) 

The decomposition is first truncated at a fixed order m and then the K m -value Gaussian vector 
constituted of the m first coordinates of the process on its Karhunen-Loeve decomposition is quantized. 
To reach optimal quantization, we have both to determine the optimal rank of truncation d x (N) (the 
quantization dimension) and to determine the optimal d x (A^)-dimensional Gaussian quantizer corre- 

d x (N) 

sponding to the first coordinates. Af(0,\f). Usual examples of such processes are the standard 

3 = 1 

Brownian motion on [0,T], the standard Brownian bridge on [0,T], the fractional Brownian motion and 
the fractional Ornstein-Uhlenbeck process. 

m 

We can also choose to use a product quantization of ^) Af(0, Xf). The product quantization is 

3=1 

the cartesian product of the optimal quantizers of the standard one-dimensional Gaussian distributions 
Af (0, ^) 1<i<d x ( N y I n the case of independent marginals, this yields a stationary quantizer. One 
advantage of this method is that the one-dimensional Gaussian quantization is a fast procedure. Newton- 
Raphson methods converge very fast to the optimal quantization (see \18\). Moreover, a sharply optimized 
database of quantizers of standard univariate and multivariate Gaussian distributions is available on the 
web site www. quantize .maths-f i . com [19| for download. Still, we have to determine quantization size on 
each dimension to obtain optimal product quantization. In this case, the minimization of the distorsion 
© comes to: 



nm,< 2^ X n™™M-Z iNn) W 2 2+ E ^n,N 1 x---xN d <N,d>l\. (4) 

n>d+l 



A solution of © is called an optimal K-L product quantizer. This problem can be solved by the "blind 
optimization procedure", which consists in computing the criterium for every possible decomposition 
N\ x • • • x Nd with N\ > ■ ■ ■ > Nd- The result of this procedure can be kept off-line for a future use. 
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Optimal decompositions for a wide range of values of N for both Brownian motion and Brownian bridge 
are available on the web site www.quantize.maths-fi.com |19| . 

In |12| . the rate of convergence to zero of the quantization error is investigated. A complete solution 
is provided for the case of Gaussian processes with regular varying eigenvalues. Rates of convergence 
are available for the above cited examples of Gaussian processes. The asymptotic of the quantization 
dimension d x (N) are investigated in |13] , The following theorem combines these results: 

Theorem 0.1 (Functional quantization asymptotics) . Let X be a centered Gaussian process on [0,T] 
with Karhunen-Loeve system (e„ , A„ ) n >i- Let (Yjv)iV>l be a sequence of quadratic optimal N— quantizers 
for X . We assume that 

A„ ~ —r as n — > oo (b > 1). 
n° 

We have: 

• span(Y/v(ft)) = spanjef , • • • ,e^ X(Ar) | and d x (N) = ft (log AT). 



. e N (X) = \\X- Y N \\ 2 ~ ^ v /6 b (&-l)-i(21ogiV)- — 
A conjecture is d x (N) ~ § log(A r ). 

It is shown in |12| that the Karhunen-Loeve eigenvalues of the fractional Brownian motion, (A^ ) n >i 
verify 

X n ~ ^2H+T as n °°> 
thus the fractional Brownian motion satifies the hypothesis of theorem 10.11 

In a constructive viewpoint, the numerical computation of the optimal quantization or the optimal prod- 
uct quantization requires a numerical evaluation of the Karhunen-Loeve eigenfunctions and eigenvalues, 
at least the very first terms. (As seen in theorem 10.11 the quantization dimension of usual Gaussian 
processes increases asymptotically as the logarithm of the size of the quantizer, so it is most likely that it 
is small. For instance, the quantization dimension d w (N) of the Brownian motion with N — 10000 is 9.) 
The Karhunen-Loeve decomposition of some usual Gaussian processes have a closed-form expression. It 
is the case of the standard Brownian motion, the Brownian bridge and the Ornstein-Uhlenbeck process. 
(The special case of the Ornstein-Uhlenbeck process is derived in [3]). 

1. The Brownian motion (Wt)te[o.T], 

er(*):=y^sin (^-1/2)1), A)T == T 1/2) ) ' , n > 1. (5) 

2. The Brownian bridge on [0,T], 

c . i . I 2 . ( t\ r. ( I 



e n (t) ■= y j, sin — J , \ n := y—j , n > 1. (6) 

3. The Ornstein-Uhlenbeck process on [0, T], starting from 0, defined by the SDE dr t = 8(mu — rt)dt + adWt, 
with a > 0, 9 > and W a standard Brownian motion on [0, T]. 



where ui\ n are the (sorted) strictly positive solutions of the equation 

6»sin(tJA„T) + ^a„ cos(cj a „T) = 0. 
4. The stationary Ornstein-Uhlenbeck process on [0, T], defined by the same SDE with ro ~ A/"(0, ao). 

e° U (t) := C n (cj\ n cos(wa„*) + 9sm(u}\ n t)) , \° U ■= — o — n>l, (8) 

where ui\ n are the (sorted) strictly positive solutions of the equation 

26»o;cos(tj A?l T) + (8 2 - uj\ n ) sin(cj All T) = 0, 
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and 

CI 2 v * v~~^" ' 2 V ' 2 ^ 7 ' 2 V" 2w a„ 

In a more general setting, we do not have a closed-form expression for the Karhunen-Loeve decom- 
position. For instance, as far as we know, the K-L expansion of the fractional Brownian motion is not 
known. Hence, a numerical method to evaluate first Karhunen-Loeve eigenfunctions is the "missing link" 
on the path to the constructive optimal quantization of more Gaussian processes. 

However, we can derive rate-optimal quantization of Gaussian processes using other series expansions 
as proposed by Luschgy and Pages in [TH[T7]. In this setting, the case of the fractional Brownian motion 
can be derived using a rate-optimal series expansion proved by Dzhaparidze and van Zanten in [7J [8]. 
Other constructive approaches for functional quantization are proposed by Wilbertz in |21j . 

In this article, we experiment the so-called "Nystrom method" [TJ [SJ [2U] for approximating the 
solution of the functional eigenvalue problem which defines the Karhunen-Loeve decomposition. First, 
we compare the result of the the numerical method with the closed-forms available for the Brownian 
motion, the Brownian bridge and the Ornstein-Uhlenbeck process. Then, the special case of the functional 
quantization of the fractional Brownian motion is handled. 

Functional quantization of Gaussian processes have numerous applications in numerical probability. 
In [3], a variance reduction method based on the functional quantization of a Gaussian process was 
proposed. This method can be seen as a "Guided Monte-Carlo simulation" (see figure [5]). Still, it 
was only applicable with Gaussian processes for which we could have a numerical evaluation of the 
Karhunen-Loeve eigenfunctions. Such a variance reduction method would be of high interest in Monte- 
Carlo simulations implying the fractional Brownian motion because its simulation schemes have a high 
complexity. 

Subsequently, we test this "functional stratification" variance reduction algorithm in option pricing 
problems within the fractional Brownian motion's counterpart of the classical Black and Scholes model. 
First, the case of a Vanilla option is benchmarked with the closed-form expression available in this case. 
Then the case of discrete barrier options is tested. 



1 The Nystrom method 

Let X be a bi-measurable Gaussian stochastic process on [0, T] defined on the probability space (f2, A, P). 
We assume that Jj Q r , E[X^]ds < oo. Let us denote Tx(t,s) the covariance function of X defined by 
Tx(t,s) — cov(X t ,X a ). The covariance operator Cx of X is defined by Cxf — Jj T , s)f(s)dt. 
It is a symmetric positive trace class operator on L 2 [0,T]. The Karhunen-Loeve basis associated with 
X, denoted (ejf) n >i is the Hilbert basis of L 2 [0,T] constituted with eigenvectors of Cx with decreasing 
eigenvalues. Now, we aim to solve numerically the eigenvalue problem 

/ Tx(-,s)f k (s)ds = X k f k , k>l. (9) 
Jo 

rp n 

The Nystrom method requires the choice of some quadrature rule L f(s)ds ~ ^ Wjf(sj). (wj)i<j< n is 

i=i 

the sequence of the weights of the quadrature rule, while (sj)i<i<n are the abscissas where / is evaluated. 
If we introduce this quadrature rule in equation ^ , we get 

n 

J2^x(t > s j )f k (s j ) = X k fk(t) te[0,T}. (10) 
Evaluating equation pop at the quadrature points yields 

n 

y^Wjrx(ti,Sj)fk(sj) = XkfkjU) i e {1, • • • ,n}. (11) 
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/ fk(tl) 

Let / be the vector 



, ((#«)) 

the matrix ((T X (U, s J ))) 1 < 4 .•<„, A = (diag(A fc )) fc= i...„ 



V fk(tn) 

and define -fT^ = KijWj. Then the eigenvalue problem becomes 

Kf = A/. (12) 

Hence, within this approximation, the functional eigenvalue problem turns into a matrix eigenvalue 
problem. As K is a covariance matrix, it is symmetric. However, since the weights are not equal for most 
quadrature rules, the matrix K is not symmetric. As outlined in |20| . numerical methods for matrix 
orthogonalization are much simpler in the symmetric case. As a consequence, we should restore the 
symmetry if possible. The method proposed in |20| is the following: 

We define the diagonal matrix D = diag(wj) and its square root D 1 ! 2 = diag( % /uJJ). Then equation 
(fl"2"|) becomes 

K-D-f = \f. (13) 

Multiplying by D 1 / 2 . we get 

[D 1/2 ■ K ■ D 1/2 ^j ■h = Xh, where h = D 1/2 ■ f. (14) 

Equation (|14p is now in the form of a symmetric eigenvalue problem. For square-integrable kernels 
(we stand in this case), this provides a good approximation of the n highest eigenvalues. 



1.1 Choice of the quadrature method 

Classical numerical methods for real symmetric matrix diagonalization are 

• The Jacobi transformation for symmetric diagonalization. 

• A tridiagonalization (by Givens or Householder reduction) followed by a QL algorithm with implicit 
shifts. 

All these numerical methods have a 0(n 3 ) complexity. As a consequence, the natural choice for the 
quadrature method would be the highest order possible (A high order Bode's formula, or a Gaussian 
quadrature) . 

However as pointed out in [11] , the Nystrom method associated with the trapezoidal integration rule 
admits an asymptotic error expansion in even powers of the step sizes as soon as the covariance function 
is differentiable (or continuous and piecewise differentiable) . As a consequence, instead of using the high 
order integration rule, we prefer to use a Richardson-Romberg extrapolation on the result of the whole 
procedure with the trapezoidal quadrature formula. We could reach an accuracy which approaches the 
machine roundoff error on the first eigenvalues when we benchmark this method on the Brownian motion, 
the Brownian bridge or the Ornstein-Uhlenbeck process. Another argument for the trapezoidal rule is 
that we encountered some small instabilities on the eigenfunction evaluation when using higher order 
schemes. 

1.2 Choice of the interpolation method 

The natural choice is to use equation (|10[) as an interpolation method for evaluating 

n 

n 

The same Richardson-Romberg extrapolation can be performed between the values of ^ WjK(t, sj)fk(sj) 

with the different orders n to compute this integral. The result is then divided by the extrapolated value 
of A fc . 

A remark on the interpolation method 
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One purpose of the quantization of a Gaussian process X, is to perform a quantization of a diffusion 
with respect X, as soon as such a stochastic integral can be defined. We can obtain a quantizer of 
the diffusion by inserting the quantizer of the Gaussian process in the diffusion equation written in 
the Stratonovich sense. The most accomplished study on this subject is [15j. In this case, we may 
also need a numerical approximation of the time-derivative of the eigenfunction in the Karhunen-Loeve 
decomposition. This work is mostly specific to the Brownian motion but main results remain valid for 
continuous semi-martingales that satisfy the Kolmogorov criteria as the Brownian bridge and Ornstein- 
Uhlenbeck processes. 

Still, a future work could be to extend these results to diffusions with respect to the fractional 
Brownian motion and other related processes. If Tx is (weakly) differentiable, a natural evaluation 

n 

method for the derivative would be f' k (t) = j- WjdiTx{t, Sj)fk(sj)- 

" 3 = 1 

One problem is that this method yields an irregular derivative. For example, this yields a piecewise 
constant derivative in the case of the Brownian motion. This causes instabilities problems when using 
Runge-Kutta integration methods for ordinary differential equations, which rely on the regularity of the 
considered Cauchy problem. 

As a consequence, a more regular interpolation method can give more satisfactory results when dealing 
with diffusions. (Spline or rational interpolation methods for instance.) 



2 Benchmark on known Karhunen-Loeve expansions 

In this section, we compare the numerical results obtained with the Nystrom methods in cases where 
we have closed-form expression of the Karhunen-Loeve expansion. The multi-steps Richardson- Romberg 
extrapolation consists in using the asymptotic error estimate of the method 

V = u n + -% + I -% + ... + o(±- 

Writing this expression for p different values of n allows us to solve a p x p linear system to nullify the 
p — 1 first orders of convergence. The three-steps Richardson-Romberg extrapolation with n = p, n = I 
and n = k gives the following solution : 

U k k 4 (m 2 - I 2 ) + Uil 4 {k 2 - to 2 ) + U m m 4 (l 2 - fc 2 ) 
(to 2 - l 2 ){l 2 m 2 + k 4 - m 2 k 2 - l 2 k 2 ) 

This result is naturally invariant by any permutation of the coefficients (k,m,l). We experienced less 
accurate results when using higher order Richardson-Romberg extrapolation, so we will settle for a 
three-steps extrapolation. 



2.1 Eigenvalues accuracy 

In tables [T] and [3J Karhunen-Loeve eigenvalues of the Brownian motion and of the Brownian bridge on 
[0, 1] are reported. Table [3] deals with the stationary Ornstein-Uhlenbeck on [0, 1] defined by the SDE 

dr t = -9r t dt + adW t , r ~ Af (o,~^ . (16) 

First column gives the theoretical value given by the closed-form. Following columns give the value 
computed with the Nystrom method with a regular step size with 25, 50 and 100 points. Last column 
gives the absolute error of a 3 steps Richardson-Romberg extrapolation method between n = 25, n = 50 
and n — 100. 

With regard to the above numerical results, Nystrom method yields a satisfactory accuracy for 
performing functional quantization of these processes. 



2.2 Eigenfunctions accuracy 

We now compare the closed-form expression of the eigenfunction with the approximation obtained by 
"Richardson-Romberg extrapolated trapezoidal Nystrom method". In tabled we report the highest 
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Figure 1: Record of the 5 highest eigenvalues of the Karhunen-Loeve decomposition of the Brownian 
motion. 
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Figure 2: Record of the 5 highest eigenvalues of the Karhunen-Loeve decomposition of the Brownian 
bridge. 
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Figure 3: Record of the 5 highest eigenvalues of the Karhunen-Loeve decomposition of the stationary 
Ornstein-Uhlenbeck process defined by the SDE dr t = —9r t dt + adWt, r$ ~ M (0, |). 



absolute difference between the closed-form expression and the approximation on a 300 points regular 
mesh of [0,1]. The tested cases are the Brownian motion, the Brownian bridge and the stationary 
Ornstein-Uhlenbeck process defined by the SDE (TIT)]) with a = 1 and 8 = 1. 

3 Quantization of the fractional Brownian motion 

The normalized fractional Brownian motion B , is a centered Gaussian process on [0, T], which has the 
following covariance function: 

r BH (t, S ) = 1 -(\tr+\ s r-\ s -tr), m 

where H e (0, 1) is called the Hurst parameter. If H = | then the process is the standard Brownian 
motion. 

A simple application of the Nystrom method presented in section [T] produces regularly shaped func- 
tional quantizers of the fractional Brownian motion. In figure a (5 x 2 x 2)— product quantizer of the 
fractional Brownian motion with 3 different values of the Hurst parameter is plotted. 
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Figure 4: Record of the biggest absolute error on the Karhunen-Loeve eigenfunctions approximation by 
the Richardson-Romberg extrapolated trapezoidal Nystrom method. The number of time steps used for 
the 3 steps interpolation are 50, 100 and 200. 300 equally spaced points on [0, 1] were tested. Each 
column corresponds to one eigenfunction. 




Figure 5: (5 x 2 x 2)— product quantizer of fractional Brownian motions on [0, 1] with Hurst exponent 
H = 0.3 (left), H = 0.5 (middle) and H = 0.7 (right). 

Still, for H < |, the covariance function of the fractional Brownian motion has singularities that 
break the convergence of the trapezoidal integration rule in even powers of the step sizes. Indeed, the 
derivative of t — > r B «(i, s) has an infinite limit for t — > + and for (t — > s~ or t — >• s + ). It breaks also 
the convergence of the whole associated Nystrom method in even powers of the step sizes. In [TJ [SI [3U], 
methods to handle such boundary and diagonal singularities are proposed. We will deal with this in 
section 13.11 

However, it is not the case for H > |, so that we can be confident in the results of this method in this 
case. In table [SI we report the 5 highest Karhunen-Loeve eigenvalues of the fractional Brownian motion 
on [0, 1] with Hurst exponent H — 0.7. The number of time steps are 128, 256 and 512. Last column 
yields the corresponding three-steps Richardson-Romberg extrapolation. All the computation has been 
performed with an octuple precision floating point number implementation to increase the accuracy of 
the 513 x 513-matrix eigensystem computation. (Let us precise that in the case of the Brownian motion 
on [0, 1], when performing the same computation, we get an absolute error smaller than le— 15 for the 
five first eigenvalues.) 

3.1 Kernel singularities when H < | 

As pointed out above, the covariance function of the fractional Brownian has a boundary singularity for 
t — » + and a diagonal singularity. In this section, we will use classical methods to handle this kind of 
singularities. See [TJ[51[5U] f° r a review of these method. 
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Figure 6: Record of the 5 highest eigenvalues of the fractional Brownian motion on [0, 1] with Hurst 
exponent H — 0.7. 



3.1.1 Handling the boundary singularity 
Change of variable 

The singular behavior of the fractional Brownian motion's covariance function V b h defined in equation 
(TIT)) can be removed by a change of variable. The change of variable u = t 2H and v = s 2H in integral 
© yields: 







h (u™,v™^j f k (v™^ —v ™ 1 dv = A fc / fc (18) 



(The second change of variable is done to preserve the symmetry of the Kernel.) 
This comes to 



rp2H 

2 



\ (M + M - I" 2 " -v^\ 2H ^j f k ^wj ~^ v ™ 1 dv = Xkfk (u^y (19) 



Quadrature rule on a single interval 

We now derive a quadrature rule on [0,T] with respect to the weight function w(v) = -^jv^" 1 = 
w(v) = 27/W Q with a :— -Jtj — 1. The aim is to make the quadrature rule exact with afnne functions as 
the trapezoidal quadrature rule is, in the case of an integration with a constant weight. 

r 1 

— x a (ax + b)dx = wi(al + b) + w r (ar + b) V(a, &) G R 2 . 
I 2i/ 

This yields 

4? (^( ra+2 ~ r+2 ) + -^r( rQ+1 " r +1 )) = a M + w - r ) + b ^ + w r) v ( fl > b ) e k 2 - 

in \a + 2 a + 1 / 

i.e. 



1 1 J \ w r 
The solution of the linear system is 

1 (a + l)l a+2 + r a+2 - (a + 2)l a+1 r 1 (a + l)r a+2 + l a+2 - (a + 2)r a+1 l 

Wl ~2H (a + l)(a + 2)(r-Z) ' Wr ~ 2H (a + \){a + 2){r - I) ' 

This is 

2277+ 1 + 2Hr2TT +1 - (2H + 1)1™ r r^r +1 + 2iJ/w+ 1 - (2H + l)r^l 



wi = 



(2H+l)(r-l) ' r (2JT+l)(r-0 



Quadrature rule for equally spaced abscissas 

Let us now consider the equally spaced abscissas points Xi — i~, i = 0, 1, • • • ,n. We now use these 
weights n times to integrate on intervals (xq , x 2H ), (x 2H , x?, H ), • • • , (x 2 J F [ 1 , ) to obtain the extended 
rule of quadrature. The convergence rate of this method is the same as the trapezoidal rule. 
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3.1.2 Handling the diagonal singularity 

We now have to handle the diagonal singularity \u — v\ 2H in equation @. One classical method if to 
use the smoothness of the solution by subtracting of the singularity. 



T BH (t,s) (f(s)-f(t))ds + r(t)f(t), 



where r(t) — J Q T B H(t, s)ds. The discretized eigenvalue problem is now transformed to 



Xkfk(U) = £ WjKij (f k {t 3 ) - f k (U)) + r(ti)f k {U) 



E WjKijfk{tj] 

3=1 



K*i) - WjKij f k (U 



(20) 



We now define the diagonal matrix D = diag(wi)i<i<„ and D 1 / 2 — diag(y / ui7)i<i<n as in section [TJ 



Moreover, we denote A = diag r(^) — £ WjKij 
Equation (|2"0"|) writes 



Ki<r 



\kfk = K-Df k + Af k . 



Multiplying by yields Xh = (d* ■ K ■ Di + a) h, with h = D^f. As a consequence, we obtain 
again a symmetric matrix eigenvalue problem. In the case of the fractional Brownian motion, the 
function r(t) = T B n(t, s)ds is derived explicitly: 



r(t) = - 



1 / T 2H+1 _ U 2H+1 



2H+1 



u 2H T _ (T-U) 



2H+1 



2H+1 



3.1.3 Optimal quantization of the fractional Brownian motion 

We now use this approximation of the Karhunen-Loeve basis to perform an optimal quantization of the 
fractional Brownian motion with a 50-100-200 three-step Richardson- Romberg extrapolated Nystrom 
method. 

In figure[71 we display the quadratic optimal N~ quantizer of the fractional Brownian motion on [0, 1] 
with Hurst exponent H = 0.25 and N — 20. In this case, the quantization dimension is 3. 




Figure 7: Quadratic TV-optimal quantizer of the fractional Brownian motion on [0, 1] with Hurst's pa- 
rameter H = 0.25 and N = 20. 
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4 Functional stratification of the fractional Brownian motion 



In this section, we experiment the functional quantization based stratified sampling algorithm proposed 
in [1] with the fractional Brownian motion. 

4.1 Background on stratification 

Let E be a separable Hilbert space. The idea of stratification is to localize the Monte-Carlo simulation 
on the elements of a measurable partition of the state space of a L 2 random variable X : (f2, A) — > (E, e). 

• Let (Ai)i(z] be a finite e-measurable partition of a E. The sets Ai are called strata. Assume that 
the weights pi = ¥(X 6 Ai) are known for i £ I and strictly positive. 

• Let us define the collection of independent random variables (X^i^i with distribution C(X\X £ 
A). 

Let F : (E,e) (K,B(R)) such that E[F 2 (X)} < +oo. 

E[F(X)] = £ E[l {XieAi} F(X)] = £ ft E[F(X)|X e A 4 ] 

iS/ iei 

= Y, Pi E[F(Xi)]. 
iei 

The stratification concept comes into play now. Let M be the global budget allocated to the computation 
of E[i^(X)] and Mi = qiM the budget allocated to compute E[i r (Xj)] in each stratus. We assume that 
Qi — 1- This leads to define the (unbiased) estimator of E[F(X)]: 

iei 

1 Mi 

f(x)m--=Y,p*w^ f{x ^ (21) 

igj 1 k=l 

where (X^)i<k<Mi is a C(X\X £ AJ-distributed random sample. 
Proposition 4.1. With the same notations: 

V 3 , v (FW I M ) = ±J2f i °h> ( 22 ) 
iei 1 

where a 2 Fi = Nsx{F(X)\X £ A t ) = Var(F(X t )) Vi £ I. 

The proof can be found in [4J. Optimizing the simulation allocation to each stratus amounts to solving 
the following minimization problem: 

min y P± a 2 where Vi = \ ( qi ) ieI € R 1 , I ft = 1 I . (23) 

In [3] , Corlay and Pages pointed out theoretical aspects of quantization that lead to a strong link between 
the problem of optimal L 2 -quantization of a random variable and the variance reduction that can be 
achieved by stratification. Three types of allocation rules for the budgets (ft)igj are proposed: 

• The "sub-optimal rule" is to set 

ft=P*, iei. (24) 

The two motivations for this choice are the facts that the weights pi are known and because it 
always reduces the variance. 

• The "optimal rule" is the solution of the constrained minimization problem (f2"3")l . The Schwartz 
inequality yields 

E( 2 2 \ 1 I 2 / \ 1 / 2 

iei Ql ) Kiel ) 

=i 
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As a consequence, the solution of the minimization problem corresponds to the equality case into 
the Schwartz inequality. Hence the solution of the minimization problem is given by 



Qi = ^ ! — ,!6 I (25) 



2 



and the corresponding minimal variance is given by Vi a F,: 

The counterpart of this method is that we do not know explicitly the solution (q*)i^j. In |10j . Etore 
and Jourdain proposed an algorithm for adaptively modifying the proportion of further drawings 
in each stratum, that converges to the optimal allocation. This can be used in a general framework. 
Another practical solution would be to implement a simple prior rough estimation of the optimal 
allocation. 

The "Lipschitz optimal" rule. When the partition (Aj)igj is a Voronoi partition associated with 
an optimal quantizer of X, Corlay and Pages considered the setting 

qi = ai, i G I, (26) 



\X - E[X\X G A t 



X G Ai 



where cr, is the local inertia of the random variable X, of 
It is proved that this setting has a uniform efficiency among the class of Lipschitz continuous 
functionals. Moreover, local inertia (crj)j E j are known. This solution overcomes the "sub-optimal 
choice" in every test done in [3]. 

4.2 On the functional stratification of Gaussian processes 

Here, we assume that X is an R- valued Gaussian process on [0, T\. We are interested in the value of 
E[F(X to , X tl , ■ ■ ■ , Xt„)] where = to < t\ < ■ ■ ■ < t n = T are n + 1 dates of interest for the underlying 
process. Let us assume that \ S Cpg(A, N) is a K-L product quantizer of X. The codebook associated 
with this product quantizer is the set of the paths of the form 



n>l 

where (e^,A^) is the Karhunen-Loeve decomposition of the process X on [0, T] and xfj 1 is the i„th 
element of an optimal quantizer of size N n of the standard one-dimensional Gaussian distribution. 
We now need to be able to simulate the conditional distribution 

C(X\X G A L ) 

where Ai is the slab associated with Xi in the codebook. 

To simulate the conditional distribution C(X\X G Aj), we will: 

• First, simulate the first K-L coordinates of X. The explicit simulation algorithm is available in [?] 

• Then simulate the conditional distribution of the marginals of the Gaussian process, its first coor- 
dinates being settled. 

In this setting, the aim is to simulate the conditional distribution 



£ Xtn , ' ' ' ,Xt. 



[ X s e?ds, [ X s e$(s)ds,--- , [ X s e^(s)ds) 
Jo Jo Jo ' 



(27) 



where (X t ) te [0,T] is a L 2 R-valued Gaussian process, and (e^, A^)fc g N* is the Karhunen-Loeve system 
associated with the process X. 

Conditional simulation: In [4], two solutions are proposed for the simulation of the conditional 
distribution (|27)) . 
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• The first one is the naive Cholesky method for Gaussian vector simulation, which has a quadratic 
complexity in the number of time steps. This first simulation scheme was not competitive for 
linearly simulable processes as the Brownian motion. In the following, we will mention this method 
as the brute force method. 

• The other solution, detailed in [3] requires a prior simulation of the unconditional distribution 
of (X to ,--- ,Xt n ) and has then a linear additional cost. This algorithm will be mentioned in 
the following as the linear conditioning algorithm. For Gaussian processes which have a linear 
simulation scheme in the unconditional case (as the Ornstein-Uhlenbeck process, the Brownian 
bridge and the Brownian motion), this method is of high interest. 

4.3 The case of the fractional Brownian motion 

Possible methods for simulating the fractional Brownian motion on a schedule to < t% < ■ ■ ■ < t n are 

• the naive Cholesky method, that has quadratic complexity, 

• and the circulant matrix method which has a 0(nln(n)) complexity [6, 22]. The circulant matrix 
method is also available for the multifractional Brownian motion [22] ■ 

No exact simulation scheme with a linear complexity exists for the fractional Brownian motion. Still, 
approximate method with linear complexity exists. If we choose the Cholesky method, there is no interest 
to use the linear conditioning algorithm proposed in |3j . The brute force Cholesky method is adapted to 
this situation. 

In every other case, if the unconditional simulation method has smaller complexity, we have interest to 
use the linear conditioning algorithm which has a linear additional cost to the unconditional simulation. 

In figure [8j we plot a few paths of the conditional distribution of the fractional Brownian motion with 
Hurst's parameter H — 0.3 knowing that they belong to a given L 2 Voronoi cell. 




Figure 8: Plot of a few paths of the conditional distribution of the fractional Brownian motion with 
Hurst's parameter H = 0.3 on [0, 3], knowing that its path belong to the L 2 Voronoi cell of the highlighted 
curve in the quantizer. 



4.4 Gaussian process reconstruction 

The first numerical test of the functional stratification of the fractional Brownian motion is a method to 
validate both the eigenfunction computation by the Nystrom method and the functional stratification 
algorithm. 

Indeed, one can rebuild the considered Gaussian process from its stratification. This yields the 
following simulation algorithm: 
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• First, simulate the discrete weighted distribution of the strata index (i,Pi)i£i to select the strata. 

• Then simulate the conditional distribution C (Xf , • • • , X tn X £ A4J of the Gaussian process in 
the strata by the method described above. 

The result should be distributed according to the distribution of the underlying Gaussian process. In 
tabled we report the covariance structure E[X ti X tj ]i<ij<n estimated by a Monte-Carlo simulation when 
X is a fractional Brownian motion with Hurst's parameter H — 0.7. The tested schedule is (i— )o<i<« 
with T = 1 and n — 5. The product decomposition of the quantization is 10 x 5 x 2. 
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1.0003 



Figure 9: Theoretical (left) and estimated (right) covariance ~E[X ti X t .] of the rebuilt fractional Brownian 
motion with H — 0.7. The number of generated paths for this Monte-Carlo simulation was 1 • 10 7 . 

In every tested case, when generating table |H1 the theoretical value lies in the 95% confidence interval. 
These confidence intervals were not displayed for briefness. We obtain the same order of accuracy with 
other values of H £ (0, 1). 



4.5 Application to option pricing 

A stochastic integral with respect to the fractional Brownian motion has been introduced in [5] by Helliot 
and van der Hoek, and in [3] by Biagini, Cksendal, Sulem and Wallner. They proposed a generalization 
of the Black-Scholes model. As in the classical Black-Scholes market, two assets are available: 



A risk-free asset whose price is given by 



dS? = rS°dt 



(28) 



• and a risky asset whose price is given by 

dS t = (iStdt + aStdB?, (29) 
where r, fi and a are constants and B H is fractional Brownian motion with Hurst parameter H . 

It has been shown that this market presents no arbitrage opportunity and is complete. Moreover, the 
solution of the stochastic differential equation (|29|) is given by 



S t = So cxp (oB? +fit- \a 2 t 2H 
The following theorem, prooved in [3] deals with the price of a European call option. 



(30) 



Theorem 4.2 (Fractional Black-Scholes Formula). The price at every time t £ [0,T] of a European call 
option with strike price K and maturity T is given by 



where 



C(t, S t ) = StAfidi) - Ke' r( - T -y^N{d 2 ) 
ln(f ) +r (T-t) + 2±(T 2H -t 2H ) 



di = 



In (§) + r(T - t) - 4(T 2ff - t 2H ) 



aVT 2H - t 2H 



(31) 

(32) 
(33) 



This closed-form expression is used to benchmark our simulation scheme of the fractional Brownian 
motion. 
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4.5.1 Benchmark with a Barrier option in the fractional Black and Scholes model 



Here, we benchmark the numerical method for a path dependent option in the case of a Barrier option 
in the fractional Black and Scholes model. For the sake of simplicity, we consider a log-normal Black 
and Scholes diffusion with no drift (no interest rate and no dividend). The chosen Hurst exponent is 
H = 0.3. The numerical results are reported in table [TU1 

The results are displayed for different values of the initial spot S, the strike K, the barrier B, the 
volatility a, the maturity T and the number of equally spaced fixing dates n. 

In this table, the first column corresponds to a simple Monte-Carlo estimator. The last three columns 
correspond to a stratified sampling estimator with different simulation allocation for each strata. 

The "sub-optimal weights" column stands for the allocation budget of equation (j2~4")h The "Lip.- 
optimal weights" column stand for the "universal stratification" budget allocation of equation (l2"tJl) . 
Both these two case have explicit allocation rules. Last column, "Optimal weights" corresponds to an 
estimation of the optimal budget allocation given in expression (|25D . 



Parameters 


Simple 
Estimator 


Strat. Estimator 
sub-optimal weights 


Strat. Estimator 
Lip. -optimal weights 


Strat. Estimator 
Optimal weights 


S = 100, K = 100 
B = 125, <t = 0.3, 
T = 1.5, n = 11 


12.5947 
[12.4429, 12.7466] 
Var = 600.5711 


12.5674 
[12.4732, 12.6615] 
Var = 230.8692 


12.5566 
[12.4654, 12.6477] 
Var = 216.3442 


12.5890 
[12.5201, 12.6579] 
Var = 123.5426 


S = 100, K = 100 
B = 200, cr = 0.3, 
T = 1, n = 11 


1.3412 
[1.2677, 1.4146] 
Var = 140.5978 


1.3826 
[1.3140, 1.4511] 
Var = 122.2808 


1.3613 
[1.3002, 1.4224] 
Var = 97.1538 


1.3769 
[1.3530, 1.4009] 
Var = 14.9352 



Figure 10: Numerical results for the Up In Call option, with 100 = x5 x 2 stratas. 



We notice that the quantization based stratified sampling method reduces noticeably the variance of the 
Monte-Carlo estimator. The universal stratification allocation rule (|26p proposed in [1] overcomes the 
sub-optimal weight allocation. Moreover, the "optimal allocation" estimation yields a better variance 
reduction factor. 

References 

[1] Kendall E. Atkinson. The numerical solution of integral equation of the second kind. Cambridge 
Monographs on Applied and Computational Mathematics, 1999. 

[2] Vlad Bally, Gilles Pages, and Jacques Printems. A quantization tree method for pricing and hedging 
multidimensional American options. Mathematical Finance, 15(1):119-168, 2005. 

[3] Francesca Biagini, Bernt Cksendal, Agnes Sulem, and Naomi Wallner. An introduction to white- 
noise theory and malliavin calculus for fractional Brownian motion. Proceedings: Mathematical, 
Physical and Engineering Sciences, 460(2041) :347-372, 2004. 

[4] Sylvain Corlay and Gilles Pages. Functional quantization based stratified sampling methods. 2010. 

[5] L.M. Delves and J.L. Mohammed. Computational methods for integral equations. Cambridge Uni- 
versity Press, 1985. 

[6] CR. Dietrich and Garry Neil Newsam. Fast and exact simulation of stationary Gaussian processes 
through circulant embedding of the covariance matrix. SIAM Journal Sci. Comput., 18:1088-1107., 
1997. 

[7] Kacha Dzhaparidze and Harry van Zanten. A series expansion of fractional brownian motion. 
Probability theory and related fields, 130:39-55, 2004. 

[8] Kacha Dzhaparidze and Harry van Zanten. Optimality of an explicit series expansion of the fractional 
brownian sheet. Statistics and probability letters, 71:295-301, 2005. 

[9] Robert J. Elliott and John van der Hoek. A general fractional white noise theory and applications 
to finance. Mathematical Finance, 13(2):301-330, 2003. 



15 



[10] Pierre Etore and Benjamin Jourdain. Adaptive optimal allocation in stratified sampling methods. 
Methodology and Computing in Applied Probability, 2008. 

[11] Han Guoqiang. Asymptotic error expansion for the Nystrom method for a nonlinear volterra- 
fredholm integral equation. Journal of Computational and Applied Mathematics, 59(1) :49 - 59, 
1995. 

[12] Harald Luschgy and Gilles Pages. Functional quantization of Gaussian processes. Journal of Func- 
tional Analysis, 196(2):486-531, December 2002. 

[13] Harald Luschgy and Gilles Pages. Sharp asymptotics of the functional quantization problem for 
Gaussian processes. Annals of Probability, 32(2), October 06 2004. 

[14] Harald Luschgy and Gilles Pages. High-resolution product quantization for Gaussian processes 
under sup-norm distortion. Bernoulli, 13(3):653-671, 2007. 

[15] Gilles Pages and Afef Sellami. Convergence of multi-dimensional quantized SDE's. 22 pages. 

[16] Gilles Pages. A space quantization method for numerical integration. J. Comput. Appl. Math., 
89:1-38, 1998. 

[17] Gilles Pages and Harald Luschgy. Expansions for Gaussian processes and Parseval frames. 2010. 

[18] Gilles Pages and Jacques Printems. Optimal quadratic quantization for numerics: the Gaussian 
case. Monte Carlo Methods and Applications, 9:135-166, 2003. 

[19] Gilles Pages and Jacques Printems. http : //www. quantize .maths-fi . com, 2005. "Web site devoted 
to optimal quantization". 

[20] William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical 
recipes in C++: The Art of Scientific Computing. Cambridge University Press, February 2002. 

[21] Benedikt Wilbertz. Construction of optimal quantizers for Gaussian measures on Banach spaces. 
PhD thesis, Universitat Trier, 2008. 

[22] Andrew T. A. Wood and Grace Chan. Simulation of stationary Gaussian processes in [0,1] d. Journal 
of Comp. and Graphical Statistics, 3:409-432, 1994. 

[23] Andrew T. A. Wood and Grace Chan. Simulation of mult ifr actional Brownian motion. Proc. Comput. 
Statist, pages 233-238, 1998. 



16 



